#	We thank Vincent Arel-Bundock and Walter Mebane for sharing code
#	with their original implementation of this approach. We thank
#	Kyle Marquardt for advice with this project.
#
#   Many thanks to Juraj Medzihorsky, who provded the base code
#	for this version
#   2017-06-03
#
#	Jonathan Hanson (jkhanson@umich.edu) contributed the remainder
#	and is responsible for any errors.
#
#
#

model {
        #   Sampling the observables:
        #       r   indexes country-years
        for (r in 1:N) {           
            #   j1 indexes variables in gaussian group 
            for (j1 in 1:n_V1) { 
                y1[r, j1] ~ dnorm(mu1[r, j1], psi[j1]) 
                mu1[r, j1] <- lambda1[j1] * x[ country[r], time[r] ] + c1[j1]
            }   
            #   j2 indexes variables in poisson group 
            for (j2 in 1:n_V2) { 
                y2[r, j2] ~ dpois( exp(mu2[r, j2]) ) 
                mu2[r, j2] <- lambda2[j2] * x[ country[r], time[r] ] + c2[j2]
            } 
        }
        #   The latents
        #       gaussian observables sampling error
        #       the unif prior on sigma is not great, used just for simplicity
        for (j1 in 1:n_V1) { 
            psi[j1] <- pow(sigma[j1], -2) 
            sigma[j1] ~ dunif(0, 1e2)
        }

        #   To help identify the model, we fix one lambda at 1 and sometimes specify 
        #	positive priors.
        
        #       indicator weights: lambdas
        for (j1 in 1:3) { lambda1[j1] ~ dnorm(0,1.0E-3) T(0,) }
        for (j1 in 4:5) { lambda1[j1] ~ dnorm(0,1.0E-3) }
        lambda1[6] ~ dnorm(0,1.0E-3) T(0,) # weberian
        lambda1[7] ~ dnorm(0,1.0E-3) T(0,) # wbstat
        lambda1[8] ~ dnorm(0,1.0E-3) # v2terr
        lambda1[9] ~ dnorm(0,1.0E-3) T(0,) # pubadmin
        lambda1[10] ~ dnorm(0,1.0E-3) # fiscalcapac
        lambda1[11] ~ dnorm(0,1.0E-3) T(0,) # infcap
        for (j1 in 12:14) { lambda1[j1] ~ dnorm(0,1.0E-3) }
        lambda1[15] <- 1
        lambda1[16] ~ dnorm(0,1.0E-3) T(0,)
        for (j2 in 1:4) { lambda2[j2] ~ dnorm(0,1.0E-3) }
        lambda2[5] ~ dnorm(0,1.0E-3) T(0,)  # AdmEffic
        for (j2 in 6:n_V2) { lambda2[j2] ~ dnorm(0,1.0E-3) }


        #       indicator intercepts: cs
        for (j1 in 1:14) { c1[j1] ~ dnorm(0,1.0E-3) }
        c1[15] <- 0
        c1[16] ~ dnorm(0,1.0E-3)
        for (j2 in 1:n_V2) { c2[j2] ~ dnorm(0,1.0E-3) }
      
        #       country-year weights: x
       for (ki in 1:N) {
            x[ country[ki], time[ki] ] ~ dnorm(0,1)
    	}    
    }
    